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ABSTRACT 

The  optimal  control  problem  of  a  typical  nuclear  reactor 
power  plant,  which  is  described  by  a  ninth-order  nonlinear 
differential  equation,  having  time-varying  parameters,  is 
considered.   The  nonlinear  model  complicates  the  optimal 
controller  synthesis.   Therefore,  the  approach  of  this  work 
is  to  approximate  the  response  of  the  reactor  system  by 
that  of  a  second-order  linear  model.   The  model  parameters 
are  chosen  to  minimize  the  derivations  between  the  system 
and  model  responses  using  a  search  routine.   The  optimal 
feedback  parameters  computed  for  the  second-order  model  5s 
used  for  suboptimal  control  of  the  system.   The  model  param- 
eters are  updated  to  reflect  the  system  nonlinearities  as 
well  as  changes  in  the  system  parameters;  the  corresponding 
control  scheme  is  adaptive.   It  is  shown  that  for  the  operat- 
ing conditions  considered,  the  adaptive  controller  need  not 
be  on-line. 

Also,  investigation  of  the  effects  of  different  weighting 
factors  in  the  cost  function,  and  the  effect  of  various  control 
rod  configurations  on  the  system  response  are  presented. 
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I.   INTRODUCTION 


Since  the  advent  of  nuclear  reactor  technology,  safety 
is  generally  a  foremost  consideration  in  the  design  of  a 
nuclear  reactor  and  its  control  system.   To  attain  the 
purpose  of  safe  operation  with  some  improvement  in  perfor- 
mance, there  is  a  possible  need  to  apply  optimal  control 
theory.   A  considerable  amount  of  research  effort  has  been 
directed  toward  the  realization  of  optimal  control  systems, 
especially  for  such  big-scale  projects  as  the  operation  of 
a  nuclear  power  reactor.   It  is  reasonable  to  assume  that 
nuclear  power  stations  of  the  future  will  be  operated 
entirely  under  the  control  of  digital  computers,  which  is 
accepted  as  the  most  suitable  device  to  perform  such  a  task. 

As  mentioned  by  Sinha  and  Bereznai  [Ref.  1],  the  over- 
simplified representation  of  the  reactor  dynamics,  usually 
the  one  delayed  neutron  group  model,  is  not  satisfactory 
for  control  purposes,  and  their  responses  deviate  consider- 
ably from  that  of  the  actual  system.   Therefore,  the  attempt 
is  made  to  evaluate  the  response  using  a  more  complete 
representation  of  the  reactor  kinetics,  and  to  include  an 
adequate  description  of  the  controller  mechanisms. 

In  this  paper,  a  realistic  system  of  a  nuclear  reactor, 
which  is  a  ninth-order  nonlinear  and  has  time-varying 
parameters,  is  used.   The  measured  quantity  to  be  controlled 
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is  the  neutron  or  power  level.   Though  any  method  that 
resolves  the  optimal  control  problem  for  this  system  is 
not  known ,  it  may  be  feasible  to  get  a  controller  that 
will  approximate  the  desired  optimal  response  sufficiently 
close  for  practical  purposes.   A  method  that  realizes 
such  a  suboptimal  performance  is  described  in  this  paper 
based  on  the  work  of  Bereznai  [Ref .  2] .   It  uses  a  second- 
order  model  to  represent  the  reactor  and  the  controller 
mechanism.   The  optimized  parameters  of  the  model  are  com- 
puted to  approximate  the  response  of  the  actual  system 
using  a  search  routine.   Then,  the  feedback  parameters 
required  to  use  the  suboptimal  feedback  controller  for  the 
nuclear  reactor  are  easily  computed.   In  this  manner,  the 
control  system  is  adapted  to  compensate  for  the  nonlinear 
system  characteristics  and  for  the  changes  in  the  system 
parameters . 
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II .   SYSTEM  REPRESENTATION 

T.he  mathematical  model  for  the  plant  under  consideration 
is  developed.   To  retain  sufficient  accuracy,  the  reactor 
kinetics  represented  by  the  six  delayed  neutron  group  is 
used.   The  unity  feedback  is  used  to  study  the  closed- 
loop  response  of  the  system. 

A.   REACTOR  KINETICS 

The  kinetic  equations  of  a  chain-reacting  pile  have 
been  derived  in  the  literature  many  times  [Ref .  3] .   The 
six  delayed  neutron  group,  space  independent,  source  free 
model  is  described  by  the  follovving  equations: 

Mtl  .  6k(t)  -  B  n(t)  +  I      X  c.Ct]       (2.1) 

i=l 

dC.(t)    6. 

-^t~=   F^nCt)  -  A.  C.(t)  (2.2) 

where  i  =  1,  2,  ...,  6. 

All  of  these  equations  are  linear  in  the  so-called  state 
variables,  n  and  C  (i  -    1,  ...,  6),  and  linear  in  the  input, 
6k(t) ,  but  the  system  is  not  jointly  linear  in  state  and 
control.   Any  feedback  of  n  to  formulate  6k(t)  results  in 
a  nonlinear  system. 
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B.   TEMPERATURE  REACTIVITY 

In  most  reactors,  heat- transfer  dynamics  are  coupled 
to  neutron  kinetics  by  a  temperature  coefficient  of  reactiv- 
ity which  is  generated  by  thermal  expansion  of  the  core  and 
by  change  in  various  neutron  cross  sections.   In  other  words, 
thermal  expansion  reduces  density,  and  thereby  the  modera- 
tion ability  and  reactivity.   Also,  the  neutron  capture-to- 
fission  ratio  of  the  fuel  is  altered  by  temperature  changes, 
and  a  smaller  effect,  due  to  changes  in  absorption  and 
fission  cross  sections,  increases  leakage  and  decreases 
reactivity.   Further,  thermal-reactivity  contribution 
results  from  so-called  Doppler  broadening  associated  with 
the  resonance  region  of  the  neutron  energy  spectrum.   The 
temperature  coefficient  a     may  usually  be  assumed  to  be  a 
constant  within  the  accuracy  with  which  it  can  be  pre- 
dicted.  Hence,  the  temperature  reactivity  may  be  expressed 
as  : 


d6k       6k     a 

+  —  n(t)  .  (2.3) 


dt       Tt    Tt 


The  significance  of  the  temperature  coefficient  of  reactivity 
is  obvious  from  Figure  2.1.   Even  if  a  linear  heat-exchange 
model  is  utilized,  the  temperature-reactivity  feedback 
causes  the  reactor  dynamical  model  to  be  nonlinear.   The 
effective  temperature  is  assumed  to  be  proportional  to  the 
operating  power  level. 
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Figure  (2.1) 
Feedback  Loop  cf  Temperature-Reactivity 

C.   ABSORBER  ROD 

The  reactivity  term  6k  is  the  sum  of  the  externally 
applied  reactivity  due  to  the  absorber  rod,  and  the  change 
in  the  reactivity  due  to  temperature  variations.   The  mathe- 
matical representation  of  the  reactivity  due  to  the  absorber 
rod  is  given  in  Reference  2: 


6k   =  6kn  +  .TW(V)  dt 
c      u    o 


(2.4) 


where  W(V)  =  0.02  V(t)     |V|  <  15 
=0.3  V   >  15 

=  -0.3  V   <  -15 


(2.5) 


that  is: 
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15 


-0.3 


V(t) 


Figure  (2.2) 
Absorber  Rod  Characteristic  Curve 

W(V)  has  the  units  mk/sec,  and  represents  the  reactivity 
change  due  to  the  motion  of  the  absorber  rod.   This  formu- 
lation assumes  that  the  control  rod  is  near  the  center  of 
its  movement,  where  the  reactivity  changes  linearly  with 
distance  from  the  center  of  the  core. 

D.   DRIVE  MOTOR 

The  time  response  of  the  drive  motor  is  given  by: 


37"  -+  -e 

dt     T       T 

m    m 


(2.6) 


where    G  is  gain,  V  is  motor  velocity, 

e  is  power  error  and  t   is  time 

v  m 

constant  of  motor. 
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E.   OVERALL  SYSTEM 

The  block  diagram  of  the  nuclear  reactor  model  considered 
in  this  thesis  is  represented  in  Figure  (2.3). 


Motor 
Drive 


Absorber 
Rod 


6  k 


Six  Group 

Neutron 
Kineti  cs 


Actual 

Power 
- — » 


5k 


Temperature 
Coef f icien" 


Figure  (2.3) 
Block  Diagram  of  Nuclear  Reactor  and  Reactivity  Mechanism 

In  order  to  study  the  closed-loop  response  of  the  system, 
unity  feedback  is  used.   The  error  signal  (e)  that  expresses 
the  deviation  between  the  demanded  and  actual  power  levels 
is  amplified,  and  the  output  is  applied  to  the  absorber  rod 
drive  motor.   The  rate  of  reactivity  insertion  (6k  /sec)  is 
proportional  to  the  error  signal  until  the  motor  reaches 
its  maximum  speed.   The  actual  reactivity  6k,  which  is  a 
measure  of  the  criticality  or  multiplication  factor  of  the 
reactor,  is  the  difference  between  the  reactivity  worth   of 
the  absorber  rod  and  the  effect  of  the  temperature  change. 
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In  state  variable  form  the  system  equations  that  represent 
the  reactor  and  the  absorber  rod  mechanism  in  an  open  loop 
configuration  are  written  in  the  following: 


fdn  \ 

/ 

dt 

dC1 

dt 

dC2 

dt 

dc3 

dt 

= 

dc4 

dt 

dc5 

dt 

dc6 

dt 

d6kt 

dt 

dV 

.  dt  J 

\ 

I        Xl      X2      X3      X4      X5      X6      I        W 


-X         000000  0 


a 


A         0         0         0         0 


0  0 


0         0      -A7      0         0         0         0  0 


0         0         0      -XA       0         0         0 


0         0         0      -A         0         0 


0         0         0         0         0 


A,       0 
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—      000000 
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10000         0000-' 
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n 
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6k. 


V 


r°  1 

0 

0 

0 

+■ 

0 

0 

0 

0 

G 

kTm 

(2.7) 


i  / 


Using  matrix  notation 

x  =  Ax  +  Be  (2.8) 

and  the  system  output  is  given  by 

y  =  Crx  (2.9) 

where  CT  =  [lp. . .0]  . 

The  numerical  values  of  the  parameters  are  given  in 
Appendix  A.   Hence,  the  parameters  of  reactor  kinetics 
equations  are  typical  values,  and  the  parameters  of  the 
absorber  rod  mechanism  are  given  in  Reference  4. 
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III.   DYNAMICAL  BEHAVIOR  OF  SYSTEM 

The  equations  of  the  system  described  in  Section  II  were 
solved  on  an  IBM  360,  Model  67  digital  computer  (W.  R.  Church 
Computer  Center  at  the  Naval  Postgraduate  School)  using  the 
method  of  digital  simulation  language  (DSL) ,  which  has  been 
suitably  modified  to  take  into  account  the  nonlinear  nature 
of  the  differential  equations. 

A.  CLOSED- LOOP  RESPONSE  OF  NUCLEAR  REACTOR 

The  closed-loop  responses  of  the  nuclear  reactor  to 
step-change  in  demanded  power  level,  for  various  initial 
power  level  are  shown  in  Figure  (3.1)  to  Figure  (3.8). 
The  reference  points  such  as  the  time  to  reach  first  full 
power  (100%  FP)  ,  the  amplitude  and  the  time  of  maximum 
overshoot  are  tabulated  in  Table  (3.1).   The  nonlinear 
characteristic  of  the  system  is  clearly  shown  in  the  figures. 

B.  EFFECT  OF  CONTROL -ROD  TO  THE  RESPONSE  OF  NUCLEAR  REACTOR 
Various  control-rod  configurations  are  tried  to  examine 

the  effect  to  the  system  response.   These  are  shown  in 

Table  (3.2).   Hence,  the  number  given  to  the  control-rod 

shown  in  Table  (3.2)  are  corresponding  to  the  number  of 

Figures  (3.9)  and  (3.10). 

The  configurations  of  1,  2  and  3  have  the  same  W 

m  d  a. 

value  and  different  slopes  for  the  linear  region.   These 
system  responses  are  shown  in  Figure  (3.9).   In  this  case, 
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Table  (3.2) 
Control  Rod  Configurations 


No 


Configuration 


Equation 


w 

mK 
-15 

0t3 

/ 

15               V 

-  n  3 

u'^      a=0.02 

W  =  0.0  2V   |V|  <  15 
=0.3     V  >  15 
=  -0.3    V  <  -15 


W 
0.3 

-5 


5 
0.3 


V 


W  =  0  .  0  6V   I V I  <  5 
=0.3     V  >  5 
=  -0.3     V  <  -5 


a=0.06 


W  =  0.01V   |V|  <  30 
=0.3     V  >  30 
=  -0.3     V  <  -30 


a=0.01 


W  =  0.0  2V   |V|  <  20 
=0.4      V  >  20 
=  -0.4    V  <  -20 


a=0.0  2 


W  =  0.02V   |V|  <  10 
=0.2      V  >  10 
=  -0.2     V  <  -10 


a=0.02 
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Figure  (3.1) 
Steady  State  Response 


10 


t (sec) 


Figure  (3.2) 
Step  Response  of  Initial  Power  Level  9 5 % 
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12  3         4 

Figure  (3.3) 
Step  Response  of  Initial  Power  Level  90% 


2         4         6         8 
Figure  (3.4) 
Step  Response  of  Initial  Power  Level  80% 


sec 
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8        12 
Figure  (3.5) 
Step  Response  of  Initial  Power  Level  601 


se< 


8         12       16     t/sec 
Figure  (3.6) 
Step  Response  of  Initial  Power  Level  50% 
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Figure  (3.7) 

Step  Response  of  Initial 

Power  Level  30% 


4         8        12 
Figure  (3.8) 
Step  Response  of  Initial  Power  Level  20% 
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the  early  system  response  for  all  three  configurations  are 

identical.   On  the  other  hand,  the  configurations  of  1,  4  and 

5  have  different  W    values.   As  understandable  from 

max 

Figure  (3.10),  the  different  W    give  much  effect  to  the 
early  time  response.   So,  the  control-rod  characteristic 
becomes  the  important  factor  in  determining  the  system 
response  for  the  early  values  of  time.   In  this  paper,  the 
control-rod  characteristic  of  1,  which  has  an  application 
in  a  particular  reactor  study  [Ref .  2] ,  is  used  in  the 
following  sections.   In  addition,  a  faster  rod  response 
characteristic  is  used  in  Section  D  where  maximum  reactivity 
insertion  is  investigated. 

C.   EFFECT  OF  TEMPERATURE  COEFFICIENT 

The  temperature  coefficient  a   is  assumed  to  be  a  constant 
but  in  a  practical  reactor  it  may  vary  over  some  range. 

To  investigate  the  effect  of  some  range  of  the  tempera- 
ture coefficient,  the  following  temperature  coefficients 
are  used  and  the  responses  are  shown  in  Figure  (3.11). 

a   =  -8  x  10""* 

a  =  -3.5  x  10"3 

a  =  1  x  10~3 

a   =  10  x  10"3 

As  the  figures  indicate,  for  values  of  the  temperature 

_  3 
coefficient  about  the  nominal  value  of  -3.5  x  10   ,  these 

temperature  coefficients  do  not  have  much  influence  to 
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the  system  response.   It  works  to  increase  or  decrease 
the  overshoot. 

D.   CONSTRAINT  OF  REACTIVITY  INSERTION 

Total  reactivity  6k  which  is  the  sum  of  the  reactivity 
due  to  the  control  rod,  and  the  temperature  reactivity  is 
expressed  as : 

6k  =  6k   +  6k  (3.1) 

In  the  previous  sections,  the  magnitude  of  the  reactivity 
was  not  constrained.   As  shown  in  Section  C,  the  magnitude 
of  temperature  reactivity  is  not  very  influential  in  the 
system  response.   The  interested  point  is  the  constraint 
applied  to  the  control-rod  reactivity  rate.   The  movement 
of  the  control-rod  should  be  restricted  to  a  maximum  value 
from  physical  viewpoints. 

Now,  the  graph  of  the  control-rod  reactivity  of  the 
initial  power  level  8  0%  is  shown  in  Figure  (3.12).   Hence, 
as  described  in  Section  B  in  Chapter  III,  the  faster 
control-rod  is  chosen.   Three  different  constraints  of 
control-rod  reactivity  are  considered.   These  three  con- 
straints are  interesting  points,  because  the  first  works 
to  reduce  both  the  first  and  the  second  overshoots,  and  the 
third  constraint  is  entirely  below  the  actual  reactivity 
except  the  starting  area.   These  situations  are  shown  in 
Figure  (3. 12)  . 
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The  Figure  (3.13)  shows  the  relation  of  the  system 
response  and  the  control-rod  reactivity.   The  obvious  fact 
is  that  the  control-rod  reactivity  gives  much  effect, 
especially  stability,  to  the  system  response.   The  influence 
of  the  constraints  of  the  control-rod  reactivity  is  shown 

in  Figure  (3.14).   As  shown  in  Figure  (3.14),  the 

_  3 
constraint  of  Sk   =  3.9  x  10    is  the  most  effective  one, 

and  the  desirable  system  response  of  fast  rise  and  minimum 

overshoot  is  obtained.   So  the  suitable  constraint  on  the 

reactivity  insertion  of  the  control  rod  contributes  to  the 

stability  of  the  system  while  improving  the  system 

performance . 
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IV.   OPTIMAL  CONTROL  PROBLEM 

There  is  a  possible  need  to  investigate  the  application 
of  optimal  control  theory  to  nuclear  reactor  system  in 
order  to  improve  the  system  performance.   In  this  thesis 
the  realization  of  optimal  control  systems  will  be  investi- 
gated for  step  demand  changes  in  power.   The  procedure  is, 
however,  applied  for  any  deterministic  demand  input. 

A.   OPTIMAL  CONTROL  LAW 

The  optimal  control  theory  has  been  derived  in  many 
literatures  [Refs .  5,  6  and  7].   The  optimal  control  problem 
is  to  find  a  control  u*eU  which  causes  the  nonlinear  system 

x(t)  =  f(x(t),  u(t),  t)  (4.1) 

to  follow  a  trajectory  x*eX  that  minimizes  the  performance 
measure 

tf 
J  =  h(x(tf),  tf)  +  /   g(x(t),  u(t),  t)  dt     (4.2) 

tQ    - 

For  most  cases  of  practical  viewpoint,  such  as  the  control 
of  a  nuclear  reactor,  the  optimal  control  must  be  realized 
in  a  feedback  manner, 

y  =  f(x(t),  t)  (4.3) 


:>b 


While  the  mathematical  solution  of  this  problem  is  established, 
it  involves  the  solution  of  a  nonlinear  partial  differential 
equation. 

Optimal -controller  synthesis  is  complicated  by  the 
nonlinear  model  which  is  utilized  to  approximate  reactor 
dynamics.   For  linear  systems,  however,  it  is  possible 
to  obtain  analytical  solutions  to  such  problems.   Further, 
it  is  shown  in  this  section  that  the  synthesis  of  the  control 
which  minimizes  an  integrated  quadratic  in  state  and  control 
only  requires  a  linear  feedback  of  state  variables. 

Assume  for  convenience  that  the  dynamical  behavior  of 
the  system  is  approximated  by  the  linear  model 

x(t)  =  A(t)  x(t)  +  B(t)  lj(t)  (4.4) 

where  x(t~)  =  x„ 

The  problem  is  to  find  the  admissible  u(t),  or,  preferably, 
u(x),  which  controls  equation  (4.4)  with  respect  to  some 
reference  x*(t)  so  as  to  minimize 


tf 

J  =  /   f0(x(t),  y(t),  t)dt  (4.5) 

o 


where 


fo(x(t),  y(t),  t)  =  |  [(x*-x)T  •  Q(t)(x*-x) 


+  u*  •  R(t)  y]  (4.6) 
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and  R(t)  is  a  real  symmetric  positive-definite  matrix  and 
Q(t)  is  a  positive-semidef inite  matrix. 

Obviously,  this  problem  may  be  solved  by  the  maximum 
principle  (or  Hamilton's  equations),  or  by  dynamic  program- 
ming (or  the  Hamilton- Jacobi-Bellman  equations) .   To  use 
the  Hamilton-Jacobi-Bellman  equation,  the  Hamiltonian  is 
formed : 

H(x(t),  u(t),  |i  ,  t) 


\    1*TM  QCt)  x(t)  +  yT(t)  R(t)  y(t)] 


T 
+  |4   •  [A(t)  x(t)  +  B(t)  u(f)]  (4.7) 


A  necessary  condition  for  y(t)  to  minimize  H  is  that 


-r—  =  0 ;  thus 

9y 


§  Cx(t),  get),  m  ,  t) 


=  R(t)  p(t)  +  BT(t)  |4  =  °         ^4-8) 


Since  the  matrix 


2 
|~4  =  RCO  (4.9) 

-y 


is  positive  definite  and  H  is  a  quadratic  form  in  y.  the 
control  that  satisfies  equation  (4.8)  does  minimize  H. 
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Solving  equation  (4.8)  for  y*(t)  gives 


y*(t)  =  -R_1(t)  BT(t)  |i  (4.10) 


which  when  substituted  in  equation  (4.7)  yields 

T 
„,  -  ^    ±s*_^      3J    .    IT,-.     1  3J  D  D-l  BT  9;J 
H(x(t),  y*(t),  g-  ,  t)  =  ^  x   Q  x  -  ^  5J-  B  R   B   ^~ 

+  t^~  A  x  (4.11) 

dx   -  ~ 


the  Hamilton- Jacobi-Bellman  eauation  is 


8J(x,t)    ,   T       n  8J(x,t)T     -.   T  3J(x,t) 
+  i  x1  Q  x  -  i-  r^— = B  R  -1  B 


3t        2  ~   -5  ~    2   3x       ~~    ~    8,x 

8J(x,t)T 
+  -^ A  X  (4.12) 

d  X         ~  ~ 

From  equation  (4.5)  the  boundary  condition  is 

J  (x(t£)  ;  tf)  =  0  (4.13) 

Since  the  minimum  cost  for  the  discrete  linear  regulator 
problem  is  a  quadratic  function  of  the  state,  it  seems 
reasonable  to  guess  as  a  solution  of  the  form  [Ref.  5] 

J    (x(t),  t)  =  \   xT  (t)  K(t)  x(t)  (4.14) 

where  K(t)  is  a  real  symmetric  positive-definite  matrix 

that  is  to  be  determined.   Substituting  this  assumed  solution 

in  equation  (4.12)  yields  the  result 
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IT*      IT        IT       -IT         T 

ix   K  x  +  |  x   Qx  -  |x   KBRX  B   Kx+x   K  A  x  =  0   (4.15) 


where  generally  all  these  matrices  and  vectors  are  functions 
of  time.   Equation  (4.15)  must  be  valid  for  all  x(t). 
Using  the  matrix  property,  this  equation  becomes: 


K(t)  +  K(t)  A(t)  +  AT(t)  K(t) 


K(t)  B(t)  R"1(t)  BT(t)  K(t)  +  Q(t)  =  0     (4.16) 


and  the  boundary  condition  is 

K(tf)  =  0  (4.17) 

Once  K(t)  has  been  determined,  the  optimal  control  law  is 
given  by 

y*(t)  =  -R_1(t)  BT(t)  K(t)  x(t)  (4.18) 

If  the  system  is  completely  controllable,  it  is  sometimes 
desirable  to  let  tr  -*■  °° ,  so  that  the  optimal  control 
accurately  maintains  the  desired  terminal  state  once  it  is 
reached.   Further,  Kalman  [Ref.  8]  has  shown  that  if 


J  =  \  f      (xT(t)  Q  x(t)  +  uT(t)  R  y(t))dt  (4.19) 


where  Q  and  R  are  positive-definite,  constant,  symmetrical 
matrices  and  the  system  [equation  (4.4)]  is  time  invariant 
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(A  and  B  are  constant  matrices) ,  then 


lim  dK(t)  a 
t+°°  dt 


(4.20) 


In  this  case,  equation  (4.16)  becomes 


ka  +  atk-kbr'1b'ik  +  q=o 


(4.21) 


For  this  problem,  the  optimal  control  is 


y*(x)  =  -R"1  BT  K  x  =  -  kT  x 


(4.22) 


and  the  synthesis  requires  a  constant  linear  feedback  of  all 
the  state  variables  as  shown  in  Figure  (4.1). 


Controller 

Figure  (4.1) 
Optimal  Feedback  Controller  for  Linear  Regulator 
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B.   APPLICATION  OF  A  SECOND-ORDER  MODEL 

The  second-order  differential  equation  is  expressed  as 


dx 


S  +  al  oT  +  a0*  =  bou 
dt 


(4.23) 


This  state  equation  is  : 


xl 

• 

= 

0 

1 

f       \ 

xl 

+ 

'0  ' 

lx2  J 

l"ao 

-au 

lx2J 

k 

u 


(4.24) 


The  elements  of  K  are  determined  by  solving  equation  (4.21) 
The  corresponding  expressions  for  the  optimal  feedback 
parameters  are  obtained  by  using  equation  (4.22),  and 
developed  in  Appendix  B  with  the  result. 


■o  +  i 


2  +  Ml  u2-,2 


bT   bQ  l    0    r    0J 


(4.25) 


ai    1     ?    ^?   ? 

t-±-  +  £-  [af  +  —  b„  +  2bnkJ 

b„    bn  L  1    r    0     0  0J 


(4.26) 


where 


Q  = 


qx   o  1 


,    R  =  r  ;    K  = 


kll   k12' 


V      V 

^  12    22J 


and 


v   _  b0k12    -     b0k22 
K0   ~~r~  '  Kl  ~r~ 
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and  selection  of  the  positive  parameters  corresponds  to  the 
positive  definite  requirement  on  k. 
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V .   OPTIMUM  APPROXIMATION  OF  NUCLEAR 
REACTOR  SYSTEM  BY  A  SECOND- ORDER  MODEL 

The  nonlinear  high-order  system  described  in  Chapter  II, 
which  is  written  below 


x  =  A  x  +  B  e  (5.1) 


where    x(t-J  =  xn  , 


e  =  x* 


T 
and      Y(t~)    =  C   x  (5.2) 


may  sometimes  be  approximated  with  sufficient  accuracy 
by  a  linear  model  for  small  variations  in  state  and  control 
To  apply  the  optimal  control  law  to  the  nuclear  reactor 
system,  this  optimum  approximation  problem  is  the  most 
necessary  thing.   This  idea  has  been  developed  by  Sinha  and 
Bereznai  [Ref.  9].   But  in  this  paper,  the  different 
error  criteria  are  used. 

A.   DESCRIPTION  OF  METHOD 

Considering  a  discrete  set  of  values  of  Yrt-\    taken  over 
a  suitable  interval  of  time, 


4  3 


Y  .■  (yj,  y2,  •  ■  • ,  yj)  (s.3) 

where  y.  =  y(ti) 

This  set  represents  discrete  sample  points  of  the  response 
of  the  system,  which  is  obtained  by  solving  equation  (5.1), 
using  the  digital  computer.   These  results  are  shown  in 
Chapter  III.   The  continuous  output  response  y^.^  is  sampled 
at  sufficiently  close  intervals  of  time  lest  the  significant 
information  is  lost.   The  objective  is  to  find  another 
output  set  y*  in  optimal  fashion,  associated  with  a  second- 
order  model  described  by  the  equations. 

x*  =  A*  x*  +  B*  u  (5.4) 

y*  =  C  x*  (5.5) 

such  that,  for  the  same  input,  the  following  objective  is 
satisfied . 

A  scalar  error  performance  function  J  is  minimized 


J  =  g[w?  (y±  -  yp]  (5.6) 


which  is  some  suitable  function  of  the  errors  (y.  -  y?)  with 
a  vector  weighting  factor  W.  attached  at  each  sampling 
instant . 

The  choice  of  the  error  criterion  expressed  by 
equation  (5.6),  has  a  direct  effect  on  the  parameters  of 
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the  approximating  model.   Since  the  purpose  of  the  error 
criterion  is  to  measure  the  extent  to  which  the  response  of 
the  model  deviates  from  that  of  the  actual  system,  the  main 
problem  is  how  to  express  this  deviation  numerically.   It 
is  normal  practice  to  consider  the  norm  of  the  output  error 
and  raise  it  to  some  power  p.   This  objective  becomes  to 
minimize  a  summation  of  deviation. 


J  =  E  W.  II  y.  -  y*  II 

1    ~i   ~i 


(5.7) 


for  the  single-output  case,  where  (y.)  is  a  scalar  sequence, 
and  the  weighting  factor  W.  is  unity,  equation  (5.7)  means 
a  measure  of  the  area  between  the  curves  when  P=l,  and  the 
mean  square  error  when  P=2. 


Response 


i+ y*(t) 


'     d  =  yCt)  -  y*(t) 


Figure  (5.1) 

Reference  Response  of  the  System 

and  the  Optimal  Second-Order  Model 
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Consider  the  system  response  Yrt-\    and  the  model  response 
y*(t\    as  shown  in  Figure  (5.1).  It  is  required  to  find  a 
model  of  the  system  so  that  the  error  criterion  is  mini- 
mized.  The  following  two  error  criteria  as  the  performance 
measure  are  used  to  approximate  the  response  y  rt\    to  the 
response  y* ,  n  of  some  model. 

(1)  minimization  of 

J±    =   I    (y   -  yp2  (5.8) 

i 

(2)  minimization  of 


J2  =  s  UXi  -  y*J2  ♦  (gi-  gi)2>       (5.9) 


by  considering  the  derivative  of  y.  and  y*,  the  more  signifi 
cant  measure  is  attained  where  the  slope  is  changing  rapidly 


B.   OPTIMUM  APPROXIMATION'  USING  PATTERN  SEARCH 

The  problem  of  approximating  high  order  system  by  a 
second-order  model  in  an  optimum  manner  can  only  be  solved 
for  specific  error  criteria  by  the  presently  available 
techniques . 

Since  an  analytical  solution  to  this  problem  does  not 
appear  feasible,  various  search  techniques  may  be  considered 
While  gradient  methods  are  generally  quite  efficient  in 
finding  a  minimum,  the  necessity  of  finding  the  partial 
derivatives  of  arbitrary  error  criteria  with  respect  to  all 
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the  model  parameters,  becomes  a  great  disadvantage. 
On  the  other  hand,  direct  search  techniques  involve 
evaluating  the  effect  of  sequential  parameter  changes  in 
an  organized  manner.   The  pattern-search  technique  of  Hooke 
and  Jeeves  [Refs.  10  and  11]  was  selected  as  a  suitable 
method  for  this  task.   Appendix  C  shows  an  outline  of  the 
pattern-search  technique,  which  is  contained  within  the 
IBM  360/67  library  under  the  subroutine  name  of  DIRECT. 

C.   INITIAL  GUESS  OF  COEFFICIENTS  OF  A  SECOND-ORDER  MODEL 
An  important  problem  associated  with  search  technique 
is  the  selection  of  the  starting  parameters,  since  these 
have  considerable  influence  on  the  convergence  of  the 
process,  and  on  the  probability  of  locating  a  local 
minimum  of  a  performance  measure.   Therefore,  the  reasonable 
initial  guess  to  compute  the  parameters  using  subroutine 
DIRECT  is  needed.   The  characteristics  of  the  second-order 
linear  system  have  been  developed  in  many  manners  [Ref .  12] . 

Since  the  response  of  a  nuclear  reactor  system  is  approxi' 
mated  by  a  second-order  model,  the  characteristic  of  that 
is  applicable  to  determine  the  coefficients.   The  differen- 
tial equation  of  a  second-order  system  in  the  form  of  the 
natural  frequency,  W  ,  and  the  damping  ratio,  £ ,  is  expressed 


as 


1_4  +  2  C  W  ^r   +  W2  x  =  KU(t)  (5.10) 

dt2     *  n  3t    n 

=  KA 
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The  response  x(t)  and  its  first  derivative  are  obtained  by 
standard  methods,  and  are  shown  below. 

x(t)  -  ^_  [i  +    X    e'^V  sin  (W  Vl    -    r}   t-*]     (5.11) 

w  2  /: — ~  n 

n        1  -  C2 


x(t)  =  p  -^ S-  sin  (W   /1  -  ?Z  t)        (5.12) 

n   ^  _  c2 


Figure  (5.2)  is  a  plot  of  x(t)  and  x(t).   Since  the  input 
is  a  unit  step,  the  value  of  A  equals  to  one.   To  make  the 
steady  state  value  of  x(t)  unity, 


x(t  )  =  ^y  =   ~-  =    1  (5.13) 

W  z   W  L 
n     n 

2 
Therefore,  if  K  =  W   ,  the  steady  state  value  of  x(t) 

becomes  unity. 

The  value  of  x    may  be  found  by  a  setting  equation 

max   J  ;  &   i 

(5.12)  to  zero,  finding  t  ,  and  substituting  this  value  of 
t   into  equation  (5.11) 


_  rW  f  

*(tD}  =  °  =  r"  ~ — ==  sin  (wn  ^  "  r;2  t}         (5-14) 

n  rx    _    c2 


Then 


W  t   =    1  (5'15) 

n  p    r 2- 

1  -  C 
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Figure    (5.2) 
x(t)    and  x(t)    vs.    time 


and 


CTT 


/ 


x(t    )    =   x  =    1    +    e 

v  p'  max 


1    -    t 


(5.16) 


for  a  special  case  o£  a  second-order  system  expressed  in 
equation  (4  .  23)  , 


x.  =  x0 
l     2 


x2  =  a0  x1   -  aa  x2  +  bQ  U(t) 


As  stated  previously,  choosing  an  =  bA  =  W   ,  a,  =  2£W 

to   0    0    n    1      n 

and  U(t)  =  unit  step  input,  the  steady  state  value  of  x..  is 
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X  (t  )  =  -»    =  1 

l       aQ 


(5.18) 


The  characteristic  response  of  the  reactor  system  of  the 
initial  power  level  901  solved  by  digital  computer  is  shown 
in  Figure  (5.3) 


100 


Power 
Level 

m 


90- 


x    =  101.44 
max 


Figure  (5.3) 
Response  of  Nuclear  Reactor  System  (90-1001) 

Substituting  t   =  2.48  and  x    =  1.0144  into  equations 
&   p  max  n 

(5.15)  and  (5.16),  the  computed  natural  frequency  and 


damping  ratio  are 


Then, 


C  =  0.8035 


W   =2.128 
n 


an  =  W  *   =  4.528 
0    n 


an  =  2£W   =3.42 
1      n 
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Using  these  values  as  the  initial  guess  coefficients,  the 
computed  optimal  parameters  by  search  routine  are: 

aQ  =  1.554 

a1  =  1.51 

There  is  considerable  difference  between  the  initial  guess 
values  and  the  computed  optimal  coefficients.   In  a  practi- 
cal point  and  considering  the  nonlinear  nature  of  the  system, 
these  initial  guess  values  are  sufficiently  near  the 
optimal  coefficients.   Since  there  is  a  definite  necessity 
to  compute  initial  guess  values,  this  method  is  shown  to 
be  acceptable . 

D.   COMPUTER  PROGRAM  OF  A  PATTERN  SEARCH 

A  computer  program  has  been  written  that  uses  a  pattern 
search  subroutine  to  find  optimum  second-order  models  for 
high-order  systems.   The  program  has  the  following  summarized 
features . 

(1)  The  output  responses  of  the  nuclear  reactor  system  to 
unit  step  function  are  available  at  discrete  uniform 
intervals  of  time. 

(2)  The  parameters  of  the  model  to  be  used  for  starting 
values  are  given  by  the  user. 

(3)  The  program  assumes  a  uniform  weighting  sequence. 

(4)  The  second-order  model  is  solved  by  the  Runge-Kutta- 
Gill  fourth-order  method.   All  calculations  are  in 
double  precision. 

(5)  The  objective  functions  to  be  minimized  are 
equations  (5.8)  and  (5.9). 
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E.  MANIPULATION  OF  THE  RESULT 

Using  the  search  routine,  the  coefficients  of  the 
approximated  second-order  model  are  manipulated  with  the 
results  shown  in  Table  (5.1),  and  the  plot  of  that  is 
Figure  [5.4) . 

The  difference  in  the  model  coefficients  using  the 
performance  measure  J1  and  that  of  the  performance  measure 
J~    is  small.   This  is  clearly  shown  in  Figure  (5.5). 

The  system  responses  and  the  approximated  second-order 
model  of  the  performance  measure  J~  to  various  initial  power 
levels  are  shown  in  Figure  (5.6)  through  Figure  (5.11). 

F.  CONSIDERATION  OF  A  THIRD-ORDER  SYSTEM 

Comparing  the  system  response  with  the  second-order 
model  as  shown  in  the  previous  part,  there  is  some  appre- 
ciable difference.   But,  as  the  considered  system  has  a 
ninth-order  nonlinearity ,  to  approximate  that  by  a  second- 
order  model,  the  difficulty  exists  in  the  natural  sense. 

So  this  difference  is  considered  acceptable.   As  a  trial 
even  though  the  optimal  control  law  of  the  third-order  model 
is  not  proposed,  the  approximated  third-order  model  is  com- 
puted using  the  same  search  routine.   The  equation  of  that 
is  : 

0+  ao  0+  ai  at +  a2x  ■  a3u^  <5-19) 
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Table  (5.1) 

Coefficients  of  Approximated  Second-Order  Model 
to  Various  Initial  Power  Level 

Initial  Power  Level  Time  Interval 


NIC  =  0.95  (28  sample  points)  0.25 

JL  =  8.191  x  10"5 
aQ  =  4.399 
ax  =  2.149 

J2  =  1.399  x  10"3 
aQ  =  4.008 
aj  =  1.809 

NIf  =  0.9  (28  sample  points)           0.25 

J   =  6.60  x  10"4 

aQ  =  1.651 

a;L  =  1.614 

J2  =  5.264  x  10~3 
aQ  =  1.554 
a1    =    1.509 

NIC  =  0.8  (28  sample  points)  0.25 

J,  =  3.792  x  10"3 
aQ  =  0.507 

a,  =  0.950 

J2  =  1.539  x  10"2 
aQ  =  0.478 
a1    =    0.885 
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N 


IC 


=  0.6 


Jn  =  1.791  x  10 


(28  sample  points) 
2 


0.5 


aQ  =  0.14 
a,  =  0.511 


J„  =  3.529  x  10 


-2 


aQ  =  0.137 
a1    =  0.496 


N 


IC 


=  0.5 


J.  =  3.705  x  10 


(28  sample  points) 
-2 


0.5 


aQ  =  0.0804 
a1  =  0.358 


J,  =  6.463  x  10 


aQ  =  0.0795 


a±    =    0.351 


N 


IC 


=  0.3 


J-,  =  7.424  x  10 


(28  sample  points) 
2 


aQ  =  2.135  x  10 
a1  =  5.689  x  10 


-2 


0.5 


1.409  x  10 


aQ  =  2.282  x  10 
a1  =  7.353  x  10 


-2 
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N 


IC 


=  0.2 


J,  =  6.097  x  10 


(28  sample  points) 
2 


0.5 


aQ  =  1.564  x  10 
a1  =  6.097  x  10 


J„  =  3.560  x  10 


-1 


aQ  =  1.633  x  10 
a,  =  6.966  x  10 
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Figure  (5.4) 
Coefficients  Plot  of  Second  Order  Model 
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Comparison  of  Step  Response  of  System 
and  Model  to  95°o  Initial  Power  Level 


Model 


100 


sec 


Figure  (5.7) 

Comparison  of  Step  Response  of  System 
and  Model  to  Initial  Power  Level  80% 
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Figure  (5.8) 

Comparison  of  Step  Response  of  System 
and  Model  to  Initial  Power  Level  60% 


sec 
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lodel 


Figure  (5.9) 

Comparison  of  Step  Response  of  System 
and  Model  to  Initial  Power  Level  501 


sec 


59 


20  t 


Figure  (5.10) 

Response  to  Initial  Power  301 
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Figure  (5.11) 
Comparison  of  Step  Response  of  System 
and  Model  to  Initial  Power  Level  20% 
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where  U(t)  =  unit  step  input.   The  obtained  data  of  the 
optimal  third-order  system  to  initial  power  level  50%  is: 

Jj  =  6.213  x  10~3 

afi  =  2 . 336 

ax  =  8.10 

a2  =  10.5 

a.  =  10.64 
J2  =  7.7  x  10"1 

aQ  =  2.16 

a  =  8.53 

a2  =  10.63 

a3  =  10.88 

The  clarified  plot  is  shown  in  Figure  (5.12).   Also,  the 
second-order  model  is  added  in  this  graph.   It  is  observed 
that  the  third-order  model  is  more  realistic  than  the 
second-order  model  to  approximate  the  actual  system.   Since 
it  is  desirable  tc  construct  the  optimal  control  law  of 
the  third-order  model,  that  work  will  be  expected  in  the 
future . 


61 


O 
Wi 


> 

o 

►J 


> 

o 

Ph 

rH 

•H 
■M 

•H 


m 

o 

0 
u) 
PJ 
O 

Pu 

/-^  to 

CNi    (1) 

rH   &< 

^— '  CD 
+■> 

a)  w 

^     >s 
3  00 

to 

•H    O 

P-4      +J 


o 
o 


0 

13 


T3 
•H 

<-> 

ci 
■H 

c 

P< 
< 


62 


VI.   SUBOPTIMAL  CONTROL  USING 


A  SECOND  ORDER  MODEL 


The  suboptimal  control  of  a  reactor  system  using  optimal 
second-order  model  is  made  possible  by  the  fact  that  the 
state  variables  available  for  feedback  are  the  actual  output 
power  level  and  its  rate  of  change.   The  optimum  linear 
feedback  of  these  variables  for  the  case  of  the  second- 
order  model  with  an  integral  quadratic  cost  function  exists, 
and  the  appropriate  feedback  coefficients  can  be  determined 
explicitly  in  terms  of  the  parameters  an  ,  a,  and  the 
constants  used  in  the  cost  function.   These  realize  a  sub- 
optimal  control  for  the  reactor  system.   The  block  diagram 
representation  of  the  suboptimal  controller,  the  reactor 
system  and  the  optimal  controller  for  the  second-order  model 
are  shown  in  Figures  (6.1)  and  (6.2). 

Using  Table  (5.1),  the  calculated  kn  and  k,  are  shown 
in  Table  (6.1).   Hence,  the  weighting  factor  q.  -  q~  =  r  in 
equations  (4.25)  and  (4.26)  is  used  to  compute  k~  and  k, . 

In  this  system  and  model  with  the  suboptimal  controller 
and  the  optimal  controller,  one  modification  is  needed 
to  make  the  steady  state  value  unity. 

From  Figure  (6.2) 


x,  (s)  a 

1 ,    ,    =  -* y (6.1) 

uTs)  2                         ^    * 

y    }  s   +  a, s  +  a0 
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Table  (6.1) 
Controller  Parameters 
to  Various  Initial  Power  Level 


Initial  Power 

Level 

ko 

95% 

0.4142 

90% 

0.4142 

80% 

0.4142 

60% 

0.4142 

50% 

0.4142 

h. 

0, 

.736 

0, 

.602 

0, 

.631 

0. 

.869 

1, 

.144 

and 


U(s)  =  k  U^s)  -  (kQ  +  kxs)  Xl(s)  (6.2) 


Substituting  (6.2)  into  (6.1) 


xx(s)  aQk 


yl(S)    s2  +  (ax  +  kx)  s  +  aQ  (1  +  kQ) 


Using  the  final  value  theorem 


(6.3) 


lim      f    ■> 
Xl(oo)     s-*0  S  XllSJ 

aQk 


aQ  (1  +  k0T  (6.4) 


Therefore,  the  feed  forward  controller  of  the  amplitude 
k  =  1  +  kn  to  make  the  steady-state  value  unity  is  needed. 
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To  make  the  comparison  possible,  between  the  responses 
of  the  suboptimal  controller  and  the  corresponding  optimal 
model  controller,  and  that  of  the  system  and  the  model 
without  any  controller,  Figure  (6.3)  shows  a  plot  of  the 
system  response  to  initial  power  level  901  without  the 
controller  and  with  the  controller.   The  describable  fact 
is  that  the  maximum  overshoot  of  the  system  response  with 
the  suboptimal  controller  is  reduced.   Also,  Figure  (6.4) 
shows  a  plot  of  the  model  response  to  the  same  initial  power 
level  without  the  controller  and  with  the  controller.   The 
desirable  reduction  of  the  time  to  reach  the  100%  power 
level  is  observed  in  the  model  response  with  the  optimal 
controller.   So,  the  considered  suboptimal  controller 
works  to  make  the  reduction  of  the  maximum  overshoot  of  the 
system  response  with  no  effect  on  the  early  time  response. 
On  the  other  hand,  the  optimal  controller  of  the  model 
does  affect  the  early  tine  response  by  making  it  faster. 
Table  (6.2)  is  the  comparison  of  the  overshoot  and  the  time 
to  reach  the  first  full  power  (1001)  of  the  system  and  the 
model,  without  controller  and  with  controller. 

The  responses  between  the  system  and  the  model  to 
various  initial  power  level  with  the  controllers  as  shown 
in  Figure  (6.5)  through  Figure  (6.9).   Comparing  these 
with  Figure  (5.5)  through  Figure  (5.9),  the  controllers  work 
to  effect  the  model  rise  time  in  such  a  manner  to  make  the 
overshoot  of  the  system  and  the  model  occur  at  the  same  time, 
and  have  approximately  the  same  value. 
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(6.3) 

Comparison  of  System  Response  of  Initial 
Power  Level  9 0%  without  Controller 
and  with  Controller 


Without 
Controller 


Figure  (6.4) 
Comparison  of  Model  Response  of  Initial 
Power  Level  90",  without  Controller 
and  with  Controller 
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Figure  (6.5) 

Step  Responses  of  System  and  Model 
with  Controllers  to  Initial  Power  Level  95% 
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Figure  (6.6) 
Step  Response  of  System  and  Model 
with  Controllers  to  Initial  Power  Level 


90% 
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Figure  (6.7) 

Step  Responses  of  System  and  Model 
with  Controllers  to  Initial  Power  Level  80% 


sec 


Figure  (6.8) 

Step  Responses  of  System  and  Model 
with  Controllers  to  Initial  Power  Level  601 
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VII.   REALIZATION  OF  ADAPTIVE  CONTROLLER 

The  effect  of  the  suboptimal  reactor  controller  and  the 

optimal  model  controller  is  described  in  the  previous 

chapter.   There  the  controller  parameters  are  fixed.   It 

is  possible  to  design    programmed  time  variations  of 

controller  parameters  to  achieve  instantaneous  near  optimum 

control  for  the  system  at  all  times.   In  addition,  the 

system  parameter  such  as  the  temperature  coefficient 

change  with  time  and  the  existence  of  such  conditions  is 

the  reason  for  application  of  an  adaptive  control  system. 

Eveleigh  [Ref.  13]  has  defined  the  adaptive  control  as 

follows : 

"If  an  index  of  performance  (IP)  is  available  which 
indicates  the  system's  instantaneous  or  short-term 
average  performance  quality,  and  if  a  control  loop 
is  set  up  to  optimize  the  IP  automatically  by  adjusting 
controller  parameters,  the  parameter-adjustment  configura- 
tion is  called  an  adaptive  control  loop.   It  is  important 
to  note  that  an  adaptive  controller  is  a  parameter- 
adjustment  loop  above  and  beyond  the  normal  feedback 
used  to  control  position,  velocity,  and  the  like. 
Adaptive  control  is  thus  an  effort  to  extend  basic 
optimum-control  concepts  to  time-varying  systems." 

A.   PROCEDURE  OF  ADAPTIVE  CONTROL 

Hence,  the  concept  of  the  adaptive  control,  which  is 
shown  in  Figure  (7.1)  as  used  in  Reference  2,  is  employed. 
This  procedure  is  described  below. 

The  model  parameters  are  precomputed,  over  the  expected 
range  of  demanded  power  level  changes,  to  be  used  as  starting 
values  for  the  optimization.   As  the  model  coefficients 
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Figure  (7.1) 
Block  Diagram  of  Adaptive  Control  System 
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are  updated,  the  parameters  of  the  feedback  controller  are 
computed  from  equations  (4.25)  and  (4.26),  depending  on  the 
nature  of  the  cost  function.   Here,  also,  the  weighting 
factor  q,  =  q-  =  r  is  used  to  make  the  comparison  possible 
in  the  previous  chapter.   That  is,  the  cost  function  has  the 
following  form: 

CO 

J  =  \  f    (x2  +  u2)  dt  (7.1) 

1    0   ~ 

Initially,  the  nuclear  reactor  is  assumed  to  be  operating 

at  a  steady  state  power  level.   At  time  t  =  0,  a  step 

change  to  a  new  operating  level  is  required.   Using  the 

model  parameters  appropriate  for  the  demanded  change  and 

the  desired  cost  function,  the  controller  parameters  are 

evaluated,  and  the  system  begins  the  transition  to  the  new 

power  level.   Over  several  sampling  intervals,  the  value 

of  the  system  response  at  each  sample  point  is  stored. 

This  means  neither  updating  of  the  model  nor  recomputation 

of  the  controller  parameters  takes  place. 

The  fitness  between  system  and  model  depends  on  the 

accuracy  of  the  starting  model.   The  model  is  reoptimized  on 

the  basis  of  the  observed  system  response.   The  next  step 

is  to  recompute  the  controller  parameters.   After  completing 

these  procedures,  the  most  optimized  controller  parameters 

of  some  interval  is  found.   This  optimized  controller 

introduces  the  desirable  system  response.   Repeating  these 

steps,  the  system  response  of  entire  intervals  is  obtained 

in  the  most  optimal  manner. 
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As  described  in  Section  B  of  Chapter  III,  the  high 
gain  of  the  control  rod  rate  gives  the  response  time 
much  influence.   To  make  faster  response,  the  following 
control  rod  is  used  instead  of  equation  (2.5). 


Figure  (7 . 2) 
High  Gain  Control  Rod  Characteristic 

The  system  response  without  controller  of  initial 
power  level  50%  using  the  high  gain  control  rod  is  shown 
in  Figure  (7.3).   Also,  the  approximated  second  order 
model  is  shown  in  the  same  plot.   The  time  to  reach  100% 
value  of  the  system  is  earlier  with  the  higher  gain  rod. 

Apparently,  judging  from  this  graph,  the  second-order 
model  does  not  fit  to  the  system  response.   Though  the 
third-order  model  describes  the  system  more  accurately 
than  the  second-order  model  as  shown  in  Figure  (5.12), 
the  usage  of  the  third-order  model  is  not  established 
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7  7 


theoretically,  and  therefore  is  not  applied  to  this 
problem. 

Figure  (7.4)  shows  a  plot  of  the  system  response  with 
the  suboptimal  controller  to  initial  power  level  501  and 
that  without  controller  to  the  same  initial  power  level. 
The  suboptimal  controller  based  on  the  optimal  control 
law  of  the  second-order  model  has  the  great  effect  to 
reduce  the  overshoot  and  the  system  response  reaches 
steady-state  value  quickly.   Apparently,  the  usage  of  the 
linear  feedback  controller  is  recommended  to  control  the 
nuclear  reactor  to  reduce  the  large  overshoot  and  to 
reach  the  steady  state  faster.   Figure  (7.5)  is  a  graph 
of  the  responses  of  the  system  and  the  model  with  the 
controllers.   Comparing  this  with  Figure  (7.3),  the 
second-order  model  with  optimal  controller  more  closely 
represents  the  reactor  system  than  when  no  controller  is 
used. 

B.   APPLICATION  OF  ADAPTIVE  CONTROL 

As  described  in  the  previous  section,  the  procedures 
of  the  adaptive  controller  are  applied  by  using  an  obser- 
vation interval  of  15  seconds  duration  and  the  adaption 
intervals  of  0.25  second.   The  obtained  parameters 
that  are  computed  repeatedly  until  sufficient  accuracy  is 
attained,  are  shown  in  Table  (7.1). 

The  corresponding  responses  of  system  and  model  with 
the  adaptive  controller  are  shown  in  Figure  (7.6).   Also, 
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Figure  (7.7)  shows  the  system  responses  with  the  non- 
adaptive  controller  and  with  the  adaptive  controller.   A 
comparison  of  Figure  (7.5)  and  Figure  (7.6)  shows  that 
there  is  some  improvement  between  the  model  response 
with  the  optimal  controller  and  that  with  the  adaptive 
controller.   In  Figure  (7.7),  there  is  no  appreciable 
difference  between  the  system  response  with  the  non- 
adaptive  controller  and  that  with  the  adaptive  controller. 
This  means  the  concept  of  adaptive  controller  is  not 
effective  as  long  as  the  system  parameters  are  constant 
and  unity  weighting  factors  are  used  in  the  cost  function. 

As  a  reference,  the  response  of  the  second-order  model 
used  to  develop  the  adaptive  controller  is  shown  in 
Figure  (7.8)  with  the  system  response  using  the  adaptive 
controller. 

To  show  the  effect  of  adaptive  controller 
clearly,  another  example  is  studied.   The  theory  of 
adaptive  control  is  applied  to  the  system  of  initial  power 
level  80%  shown  in  Figure  (7.9).   The  tabulated  result 
is  shown  in  Table  (7.2). 

Figure  (7.9)  is  a  plot  of  the  system  responses  with 
the  non-adaptive  controller,  and  with  the  adaptive  con- 
troller.  Though  there  is  very  little  difference  between 
them,  they  are  almost  considered  identical. 
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VIII.   INFLUENCE  OF  COST  FUNCTION 


As  having  been  described,  the  cost  function  is: 


00 

J  =  \  S    (xT  Q  x  +  ry2)  dt  (8.1) 

z  0   ~ 


In  the  previous  chapter,  the  controller  parameters  were 
computed  for  the  case  that  Q  is  an  identity  matrix  and  r 
equals  to  unity. 

In  this  chapter,  the  various  weighting  factors  for 
the  cost  functions  are  considered,  depending  on  the  optimal 
control  law.   The  selection  of  appropriate  values  for  Q 
and  r  is  important,  and  goes  to  the  core  of  the  optimal 
regulator  problem.   The  process  of  selection  consists  of 
assigning  certain  values  to  Q  and  r,  incorporating  the 
resulting  controller  into  the  system. 

Equations  (4.25  and  (4.26)  of  the  controller  parameters 

are  rewritten  for  convenience 

1 
i       a0    1   r  2    ql  ,2,2 

0  Dq   b„  L  0   r   0J 

1 

K    =  "  xr-   +  J-  [a?  +  —  bn  +  2bnkJ2 

1  bT   bQ  L  1   r   0     0  0J 

where  q1  and  q?  are  the  components  of  a  matrix  q,  r  is  the 
relative  weighting  factor. 

These  weighting  factors  of  the  cost  function  decide 
the  controller  parameters. 
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The  following  values  are  tried  to  check  the  behavior 
of  the  actual  system. 


ql   q2 


The  computed  controller  parameters  are 


ql   q2 
(2)  —  =  —  =  1 

v  *    r    r 


ql   q2 
(3)  —  =  —  =  10 
K    J    x  r 


ql    q2 
(4)  _k  =  _A  =  20 

v    r    r 


kQ  =  0.0488 
k±   =  0.1426 


kQ  =  0.4142 
k-  =  0.797 


kQ  =  2.317 
kx  =  3.008 


kQ  =  3.583 
kx  =  4.558 


These  responses  are  shown  in  Figure  (8.1).   The  stability 
of  the  system  with  the  suboptimal  controller  depends  on 
weighting  factors  of  the  cost  function  apparently.   Since 
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the  optimal  control  law  is  based  on  the  quadratic  cost 
function,  it  is  essential  to  select  the  most  reasonable 
cost  function,  as  described  in  Reference  14.   The 
resulting  fact  is  consistent  with  Reference  14. 
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IX.   CONCLUSIOXS  AND  RECOMMENDATIONS 

The  following  conclusions  and  recommendations  are 
offered  as  a  result  of  the  investigation  contained  in 
this  report. 

A.   CONCLUSIONS 

1.  The  simulation  of  the  dynamic  behavior  of  the 
nuclear  reactor  system  on  the  digital  computer  to  make 
feasible  this  investigation  was  obtained  with  some 
difficulty. 

To  obtain  the  stabilized  system  response  required 
considerable  effort,  because  the  system  is  nonlinear  and 
described  by  a  ninth-order  differential  equation.   The 
unknown  factors  such  as  the  specified  values  of  the 
nuclear  reactor  parameters  complicated  the  difficulty  of 
implementation  of  the  system  equations. 

2.  In  this  paper,  the  gain  of  the  control  rod,  the 
reactivity-rate  constraint  and  the  effect  of  the  tempera- 
ture coefficient  have  been  studied.   The  high  gain  of  the 
control  rod  causes  the  faster  response,  and  the  reactivity- 
rate  constraint  suitably  chosen  makes  possible  the  more 
desirable  response  regarding  the  stability,  particularly 
to  reduce  the  deviation  from  the  desired  power  level.   The 
magnitude  of  the  temperature  coefficient  has  a  lesser 
effect  on  the  system  response. 


92 


3.  It  has  been  shown  that  the  results  from  the 
optimal  second-order  model  of  the  reactor  may  be  used 
effectively  for  suboptimal  control  of  the  reactor  system. 
This  suboptimal  controller  of  the  reactor  reduces  the 
overshoots  of  the  system  response.   Considering  these 
overshoots  are  not  desirable  practically,  its  usage 
provides  an  improvement  in  system  performance.   To  find 
the  optimal  second-order  model,  the  computer  program  of 
search  technique  is  applied  to  obtain  the  second-order 
model  parameters. 

4.  An  adaptive  control  system,  which  may  be  realized 
in  practice  to  control  the  power  level  changes  of  a 
nuclear  reactor,  has  been  proposed. 

The  parameters  of  the  second-order  linear  model 
are  continuously  updated  so  that  the  model  accurately 
represents  the  behavior  of  the  system.   Also,  the  parameters 
of  the  controller  are  updated  to  achieve  the  concept  of 
the  adaptive  control.   The  result  of  this  study  shows  no 
improvement  of  the  system  response  using  the  adaptive 
controller  for  the  operating  condition  considered.   The 
adaptive  control,  therefore,  would  not  have  to  be  on- 
line, and  could  be  used  to  update  the  model  parameters 
to  account  for  slowly  changing  parameters  of  the  reactor 
system. 

5.  The  weighting  factor  of  the  cost  function  gives 
the  stability  of  the  system  response  much  influence. 
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Based  on  physical  considerations,  the  weighting  factor 
should  be  chosen  as  correctly  as  possible  to  provide  the 
desired  response  and  this  would  be  a  trial  and  error 
procedure . 

B.   RECOMMENDATIONS 

Further  workers  in  this  area  should  consider  the 
following  facets  of  nuclear  reactor  systems. 

1.  The  nuclear  reactor  system  studied  in  this  report 
did  not  consider  the  noise  problem.   For  the  purpose  of 
the  present  paper,  it  was  assumed  that  an  instantaneous 
noise-free  measure  of  the  reactor  power  level  is  available. 
However,  as  the  temperature  channel  signal  indicates  the 
power  level  of  the  reactor  at  an  earlier  instant,  because 
of  the  finite  transport  tine  between  the  reactor  core 

and  the  temperature  transducer,  this  signal  has  a  con- 
siderable noise  component  due  to  the  turbulent  coolant 
flow. 

Therefore,  the  noise  consideration  is  strongly 
recommended  for  future  studies. 

2.  Since  the  proposed  scheme  relies  on  the  use  of 
search  routines,  improvement  in  the  efficiency  of  the 
method  particularly  as  applied  to  digital  process  computers, 
would  greatly  enhance  the  usefulness  of  the  new  technique. 

3.  In  this  report,  the  effectiveness  of  the  adaptive 
control  has  not  been  shown.   Consideration  should  be 
given  to  the  adaptive  controller  for  different  weighting 
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factors  in  the  cost  function  to  find  out  whether  an  on-line 
adaptive  scheme  would  be  required. 

4.   In  much  of  optimal  control,  the  natural  formulation 
is  by  state  variables.   A  well-known  and  popular  result  is 
that  for  the  linear  regulator  problem  with  quadratic 
performance  index  as  applied  in  this  study.   It  has  been 
shown  [Ref.  15],  however,  that  the  solution  is  far  from 
realistic  or  optimum  in  an  engineering  sense.   Nevertheless, 
the  results  and  scheme  obtained  for  the  linear  quadratic- 
index  regulator  problem  have  been  forced  indiscriminately 
on  the  linear  time- invariant  problem. 

Since  the  weakness  of  the  state-variable  formula- 
tion in  coping  with  linear  feedback  control  system  design 
is  especially  apparent  in  the  sensitivity  problem,  which 
is  one  of  the  primary  reasons  for  the  use  of  feedback,  the 
sensitivity  theory  by  Horowitz  and  Shaked  [Ref.  16]  should 
be  investigated. 
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APPENDIX  A 
NUMERICAL  VALUES  OF  THE  REACTOR  MODEL  PARAMETERS 


3  =  6.41  x  10"3 

3X  =  2.70  x  10"4 

32  =  7.40  x  10"4 

33  =  2.53  x  10"3 

34  =  1.26  x  10"3 

35  =  1.40  x  10"3 

36  =  2.10  x  10~4 
A-,  =  3.014  sec 
^2  ~  1 • 136  sec 
X3  =  0.301  sec"1 

A4  =  0.111  sec"1 

-  2     - 1 

Xr  =  3.05  x  10    sec 

X6  =  1.24  x  10"2  sec"1 

I  =  10"3  sec 
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xm  =  0.16  sec 

t.  =  12.5  sec 

-3 
a   =  -3.5  x  10 
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APPENDIX  B 
DERIVATION  OF  EQUATIONS  ON  kQ  AND  k 


The  equation  is 


where 


KA+ATK-KBR_1BTK+Q=0 


A  = 


I  ~ao     ~au 


0 


B    = 


Q  ■ 


0     1 


R   =    r    ,    K   = 


11  12 


21  22 


Computing  term  by  term 


with  k12  =  kn 


r 


K  A  = 


a0k12     kll  "  alk12>l 


I  'a0k22 


k21  "  aik22J 


AT  K 


"a0k21 


"a0k22 


kll  '  alk21    k12  "  alk22 


-]   T 
B  R  J  B1 


70     0 
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-IT        0 
K  B  R  x  B1  K  =  — - 


k12k21 


k22k21 


k12k22 


22 


After  substituting  these  matrices  into  the  equation, 
rearranging,  the  following  four  equations  are  obtained 

u2 


~a0  k12  "  a0  k12 


-T-   k12  +  «1  "  ° 


-a~  k„„  +  k.. .,  -  a,  k 


0 


k,  „  k_  =  0 


0  "22     11     1   12    r   "12   22 


k.,  n  -  a.,  k01  -  an   k 


k, „  k„„  =  0 


ll     l   21     0  "22     r   "12   22 


a,  k„„  +  k^.,  -  a-,  k 


12     1   22    "21    "1   22    r   "22   M2 


0,2         n 
k99  +  q   =  0 


From  the  first  one  of  the  above  four  equations 


T2  k12  +  2a0  k12  "  «»1  "  ° 


So , 
k12 


2  +  b  0 

ao  \  ao  +  —  qi 


a0r    r   I  2   *1  .  2 

+  an  +  —  b  n 

b  0   b  o 


Qfc 


Also, 


TT  k22    +    2alk22    "    (2k12    +   ^    =    ° 


2    .    b'( 


1    +    Jal    +   —    C2k12    +    «2) 


22 


From  the  optimal  control  law, 


-IT         T 
U*  =  -  R   B   K  X  =  -k  x 


then, 


kT  =  [kn  kj  =  R"1  BT  K  =  1  [0 


0  "1 


0 


'   k      k   ^ 
Kll    K12 


V      k 

21     22  ) 


~   F~  [k21   k22] 


So , 


0     r 


,   =  b0k22 
1     r 
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Therefore 


k0    = 


■o   .    1 


[a 


1   ,2,7 


^       bQ    ^0        r      u0 


J  +  B^  [ai  +  T-  b2o  +  2boV 
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APPENDIX  C 
OUTLINE  OF  THE  PATTERN  SEARCH 

As  devised  by  Hooke  and  Jeeves  [Ref .  10] ,  pattern 

search  is  a  direct  search  method  which  results  in  relatively 

efficient  search  along  straight  ridges  or  ravines.   It  is 

attempted  with  pattern  search  to  establish  the  pattern 

of  successful  search  points  in  the  immediate  past  from 

which  plausible  future  search  points  are  predicted.   The 

method  for  the  two-dimensional  case  is  illustrated  in 

Figure  (C.l).   The  search  progresses  from  an  initial  base 

point  x° ,  and  a  minimum  of  f(x)  is  to  be  found  without  the 

use  of  derivatives  of  f(x).   A  small  displacement  d-, 

from  x°  in  the  x,  direction  is  effected,  and  f(x°  +  d,  ,  x°) 

is  evaluated  and  compared  with  f(x°   x°)  .   If  the  latter 

is  smaller  than  the  former,  a  small  displacement  in  the 

-X-.  direction  is  effected,  and  f  (x?  -  d,  ,  xp  is  compared 

with  f(x°   x°).   If  the  latter  is  greater  than  the  former, 

a  small  displacement  d?  is  made  from  (x°  -  d, ,  x°)  in  the 

x2  direction,  and  f(x°  -  d, ,  x°)  is  compared  with  f (x°  -  d, , 

x^  +  d?).   Assuming  the  latter  is  the  smaller  of  the  two, 

a  pattern  move  is  made  in  the  direction  established  by  the 

preceding  successful  moves  to  the  point  x   in  Figure  (3.1). 

At  x  ,  the  small  moves  which  proved  to  be  successful  around 

x°  are  repeated  and  if  successful  lead  to  a  second  pattern 
2  2 


mov 


e  to  x  .   From  x  ,  small  displacements  of  magnitude  d. 
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Figure  (CI) 


Pattern  Search 
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along  the  x..  coordinate  prove  to  be  unsuccessful,  whereas 

2    2  2    2 

the  displacement  d2  yields  f (x1  ,  x2  +  d-J  <  f (x^  ,  x  ) . 

A  pattern  move  is  made  therefore  along  the  line  which 

passes  through  the  points  (x,  -  d..  ,.  x?  +  d~)  and 

(x^  ,  x22   +  d2). 

Note  that  the  step  size  of  the  pattern  search  is  made 

progressively  larger  as  the  general  directions  dictated 

by  the  pattern  prove  successful.   When  failure  occurs, 

f 
as  at  point  x  ,  the  step  size  of  the  pattern  move  is 

3       3 
reduced  as  depicted,  to  point  x  .   At  x  ,  the  displacement 

3 
d,  from  x  along  the  x,  coordinate  fail  to  improve  f(x),  as 

3  2 

do  displacement  d?  from  x   along  the  x   coordinate.   Reduc- 
tions in  d-,  and  d~  are  therefore  required  at  this  point, 
and  the  process  starts  anew.   Ultimately  both  the  dis- 
placement and  the  pattern  step  sizes  are  reduced  below 
some  preassigned  limit,  and  the  search  is  terminated. 
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